import numpy as np
from pyproj import Proj, Transformer

def estimate_base_station(ksxt_data):
    lat_sum, lon_sum, count = 0, 0, 0
    utm_sum_x, utm_sum_y = 0, 0
    
    # 定义WGS84和UTM转换
    wgs84_to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32650", always_xy=True)  # UTM Zone 50N (示例)
    utm_to_wgs84 = Transformer.from_crs("EPSG:32650", "EPSG:4326", always_xy=True)
    
    for line in ksxt_data:
        parts = line.split(',')
        if len(parts) < 16:
            continue
        
        try:
            lon = float(parts[2])  # 车载GPS经度
            lat = float(parts[3])  # 车载GPS纬度
            easting_offset = float(parts[14])  # 东向偏移（米）
            northing_offset = float(parts[15])  # 北向偏移（米）
            
            # 先转换流动站的经纬度到 UTM
            utm_x, utm_y = wgs84_to_utm.transform(lon, lat)
            
            # 计算基准站的 UTM 坐标
            base_utm_x = utm_x - easting_offset
            base_utm_y = utm_y - northing_offset
            
            # 累加计算基准站的 UTM 坐标
            utm_sum_x += base_utm_x
            utm_sum_y += base_utm_y
            count += 1
        except ValueError:
            continue
    
    if count == 0:
        return None, None, None, None
    
    avg_utm_x = utm_sum_x / count
    avg_utm_y = utm_sum_y / count
    
    # 转换回 WGS84 经纬度
    base_lon, base_lat = utm_to_wgs84.transform(avg_utm_x, avg_utm_y)
    
    return base_lat, base_lon, avg_utm_x, avg_utm_y

# 示例数据（请替换为实际数据）
ksxt_lines = [
    "$KSXT,20250305012108.60,116.35054828,39.99146561,35.5831,32.38,0.54,340.30,0.012,0.00,3,3,32,34,-1.623,16.138,0.001,-0.004,0.012,-0.002,,*00",
    "$KSXT,20250305012108.65,116.35054828,39.99146562,35.5850,32.42,0.63,346.61,0.015,0.00,3,3,32,34,-1.623,16.139,0.002,-0.004,0.015,-0.002,,*08",
    "$KSXT,20250305012109.00,116.35054829,39.99146563,35.5823,32.33,0.69,339.24,0.015,0.00,3,3,32,34,-1.623,16.140,-0.000,-0.005,0.014,-0.002,,*2A"
]

base_lat, base_lon, base_utm_x, base_utm_y = estimate_base_station(ksxt_lines)
print(f"Estimated Base Station Coordinates: Latitude = {base_lat}, Longitude = {base_lon}")
print(f"UTM Coordinates: Easting = {base_utm_x}, Northing = {base_utm_y}")
